{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1: Import Networks from Statoil Files\n",
    "This example explains how to use the OpenPNM.Utilies.IO.Statoil class to import a network produced by the Maximal Ball network extraction code developed by Martin Blunt's group at Imperial College London. The code is available from him upon request, but they offer a small library of pre-extracted networks on their [website](https://www.imperial.ac.uk/engineering/departments/earth-science/research/research-groups/perm/research/pore-scale-modelling/micro-ct-images-and-networks/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:32:45.929737Z",
     "iopub.status.busy": "2021-06-24T11:32:45.928107Z",
     "iopub.status.idle": "2021-06-24T11:32:46.656232Z",
     "shell.execute_reply": "2021-06-24T11:32:46.654702Z"
    }
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "import scipy as sp\n",
    "import numpy as np\n",
    "import openpnm as op\n",
    "%config InlineBackend.figure_formats = ['svg']\n",
    "np.set_printoptions(precision=4)\n",
    "np.random.seed(10)\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following assumes that the folder containing the 'dat' files is in a directory called 'fixtures' in the same directory as this script.  You can also enter a full path to the files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:32:46.664804Z",
     "iopub.status.busy": "2021-06-24T11:32:46.661474Z",
     "iopub.status.idle": "2021-06-24T11:32:46.771710Z",
     "shell.execute_reply": "2021-06-24T11:32:46.772925Z"
    }
   },
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "path = Path('../_fixtures/ICL-Sandstone(Berea)/')\n",
    "project = op.io.Statoil.import_data(path=path, prefix='Berea')\n",
    "pn = project.network\n",
    "pn.name = 'berea'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This import class extracts all the information contained in the 'Statoil' files, such as sizes, locations and connectivity. Note that the ``io`` classes return a ``project`` object, and the network itself can be accessed using the ``network`` attribute.  The following printout display which information was contained in the file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:32:46.781101Z",
     "iopub.status.busy": "2021-06-24T11:32:46.779409Z",
     "iopub.status.idle": "2021-06-24T11:32:46.788733Z",
     "shell.execute_reply": "2021-06-24T11:32:46.789498Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "openpnm.network.GenericNetwork : berea\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "#     Properties                                    Valid Values\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "1     pore.area                                      6298 / 6298 \n",
      "2     pore.clay_volume                               6298 / 6298 \n",
      "3     pore.coords                                    6298 / 6298 \n",
      "4     pore.radius                                    6298 / 6298 \n",
      "5     pore.shape_factor                              6298 / 6298 \n",
      "6     pore.volume                                    6298 / 6298 \n",
      "7     throat.area                                   12098 / 12098\n",
      "8     throat.clay_volume                            12098 / 12098\n",
      "9     throat.conduit_lengths.pore1                  12098 / 12098\n",
      "10    throat.conduit_lengths.pore2                  12098 / 12098\n",
      "11    throat.conduit_lengths.throat                 12098 / 12098\n",
      "12    throat.conns                                  12098 / 12098\n",
      "13    throat.length                                 12098 / 12098\n",
      "14    throat.radius                                 12098 / 12098\n",
      "15    throat.shape_factor                           12098 / 12098\n",
      "16    throat.total_length                           12098 / 12098\n",
      "17    throat.volume                                 12098 / 12098\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "#     Labels                                        Assigned Locations\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "1     pore.all                                      6298      \n",
      "2     pore.inlets                                   201       \n",
      "3     pore.outlets                                  246       \n",
      "4     throat.all                                    12098     \n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n"
     ]
    }
   ],
   "source": [
    "print(pn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point, the network can be visualized in Paraview. A suitable '.vtp' file can be created with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:32:46.822698Z",
     "iopub.status.busy": "2021-06-24T11:32:46.795885Z",
     "iopub.status.idle": "2021-06-24T11:32:46.983377Z",
     "shell.execute_reply": "2021-06-24T11:32:46.982027Z"
    }
   },
   "outputs": [],
   "source": [
    "op.io.VTK.export_data(network=pn, filename='imported_statoil')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting network is shown below:\n",
    "\n",
    "<img src=\"http://i.imgur.com/771T36M.png\" style=\"width: 60%\" align=\"left\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clean up network topology\n",
    "\n",
    "Although it's not clear in the network image, there are a number of isolated and disconnected pores in the network.  These are either naturally part of the sandstone, or artifacts of the Maximal Ball algorithm.  In any event, these must be removed before proceeding since they cause problems for the matrix solver.  The easiest way to find these is to use the ```check_network_health``` function on the network object.  This will return a dictionary with several key attributes, including a list of which pores are isolated. These can then be trimmed using the ``trim`` function in the ``topotools`` module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:32:46.990482Z",
     "iopub.status.busy": "2021-06-24T11:32:46.989079Z",
     "iopub.status.idle": "2021-06-24T11:32:47.003766Z",
     "shell.execute_reply": "2021-06-24T11:32:47.004912Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of pores before trimming:  6298\n",
      "Number of pores after trimming:  6004\n"
     ]
    }
   ],
   "source": [
    "print('Number of pores before trimming: ', pn.Np)\n",
    "h = pn.check_network_health()\n",
    "op.topotools.trim(network=pn, pores=h['trim_pores'])\n",
    "print('Number of pores after trimming: ', pn.Np)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dealing with Inlet and Outlet Pores\n",
    "When importing Statoil networks, OpenPNM must perform some 'optimizations' to make the network compatible.  The main problem is that the original network contains a large number of throats connecting actual internal pores to fictitious 'reservoir' pores.  OpenPNM strips away all these throats since 'headless throats' break the graph theory representation.  OpenPNM then labels the real internal pores as either 'inlet' or 'outlet' if they were connected to one of these fictitious reservoirs.  \n",
    "\n",
    "It is fairly simple to add a new pores to each end of the domain and stitch tehm to the internal pores labelled 'inlet' and 'outlet', but this introduces a subsequent complication that the new pores don't have any geometry properties.  For this example, we will not add boundary pores, but just the pores on the inlet and outlet faces.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: Calculating Permeability of the Network\n",
    "\n",
    "### Setup Geometry, Phase, and Physics Objects\n",
    "\n",
    "In OpenPNM 2+ it is optional to define Geometry and Physics objects (These are really only necessary for simulations with diverse geometrical properties in different regions, resulting in different physical processes in each region, such as multiscale networks for instance).  It is still necessary to define **Phase** objects:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:32:47.012030Z",
     "iopub.status.busy": "2021-06-24T11:32:47.010305Z",
     "iopub.status.idle": "2021-06-24T11:32:47.097202Z",
     "shell.execute_reply": "2021-06-24T11:32:47.095933Z"
    }
   },
   "outputs": [],
   "source": [
    "water = op.phases.Water(network=pn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Apply Pore-Scale Models\n",
    "\n",
    "We must add the hagen-poiseuille model for calculating the conductance.  In OpenPNM 2+ it is possible to add Physics models to Phase objects, which is often simpler than than applying the same model to multiple Physics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:32:47.103715Z",
     "iopub.status.busy": "2021-06-24T11:32:47.102300Z",
     "iopub.status.idle": "2021-06-24T11:32:47.112108Z",
     "shell.execute_reply": "2021-06-24T11:32:47.113193Z"
    }
   },
   "outputs": [],
   "source": [
    "water.add_model(propname='throat.hydraulic_conductance',\n",
    "                model=op.models.physics.hydraulic_conductance.valvatne_blunt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that boundary pores and throats had no geometrical properties associated with them, so the   hydraulic conductances of boundary throats will be undefined (filled with NaNs):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:32:47.121089Z",
     "iopub.status.busy": "2021-06-24T11:32:47.119098Z",
     "iopub.status.idle": "2021-06-24T11:32:47.125299Z",
     "shell.execute_reply": "2021-06-24T11:32:47.124216Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.1825e-13 1.6357e-13 4.3332e-14 ... 1.6230e-11 7.9766e-12 6.7575e-11]\n"
     ]
    }
   ],
   "source": [
    "print(water['throat.hydraulic_conductance'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run StokesFlow Algorithm\n",
    "Finally, we can create a **StokesFlow** object to run some fluid flow simulations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:32:47.133935Z",
     "iopub.status.busy": "2021-06-24T11:32:47.132771Z",
     "iopub.status.idle": "2021-06-24T11:32:47.259834Z",
     "shell.execute_reply": "2021-06-24T11:32:47.259159Z"
    }
   },
   "outputs": [],
   "source": [
    "flow = op.algorithms.StokesFlow(network=pn, phase=water)\n",
    "flow.set_value_BC(pores=pn.pores('inlets'), values=200000)\n",
    "flow.set_value_BC(pores=pn.pores('outlets'), values=100000)\n",
    "flow.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting pressure field can be visualized in Paraview, giving the following:\n",
    "\n",
    "<img src=\"https://i.imgur.com/AIK6FbJ.png\" style=\"width: 60%\" align=\"left\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Determination of Permeability Coefficient\n",
    "\n",
    "The way to calculate K is the determine each of the values in Darcy's law manually and solve for K, such that $$ K = \\frac{Q\\mu L} {\\Delta P A} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:32:47.264958Z",
     "iopub.status.busy": "2021-06-24T11:32:47.264378Z",
     "iopub.status.idle": "2021-06-24T11:32:47.270645Z",
     "shell.execute_reply": "2021-06-24T11:32:47.270136Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.1566e-12]\n"
     ]
    }
   ],
   "source": [
    "# Get the average value of the fluid viscosity\n",
    "mu = np.mean(water['pore.viscosity'])\n",
    "# Specify a pressure difference (in Pa)\n",
    "delta_P = 100000\n",
    "# Using the rate method of the StokesFlow algorithm\n",
    "Q = np.absolute(flow.rate(pores=pn.pores('inlets')))\n",
    "# Because we know the inlets and outlets are at x=0 and x=X\n",
    "Lx = np.amax(pn['pore.coords'][:, 0]) - np.amin(pn['pore.coords'][:, 0])\n",
    "A = Lx*Lx  # Since the network is cubic Lx = Ly = Lz\n",
    "K = Q*mu*Lx/(delta_P*A)\n",
    "print(K)"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
